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ABSTRACT 

We examine the bar instability in galactic models with an exponential disk and a cuspy dark matter 
(DM) halo with a Navarro-Frenk- White (NFW) cosmological density profile. The equilibrium models 
are constructed from a 3-integral composite distribution function but subject to the bar instability. 
We generate a sequence of models with a range of mass resolution from 1.8K to 18M particles in the 
disk and 10K to 100M particles in the halo along with a multi-mass model with an effective resolution 
of ~ 10 10 particles. We describe how mass resolution affects the bar instability, including its linear 
growth phase, the buckling instability, pattern speed decay through the resonant transfer of angular 
momentum to the DM halo, and the possible destruction of the halo cusp. Our higher resolution 
simulations show a converging spectrum of discrete resonance interactions between the bar and DM 
halo orbits. As the pattern speed decays, orbital resonances sweep through most of the DM halo 
phase space and widely distribute angular momentum among the halo particles. The halo does not 
develop a flat density core and preserves the cusp, except in the region dominated by gravitational 
softening. The formation of the bar increases the central stellar density and the DM is compressed 
adiabatically increasing the halo central density by 1.7 x. Overall, the evolution of the bar displays a 
convergent behavior for halo particle numbers between 1M and 10M particles, when comparing bar 
growth, pattern speed evolution, the DM halo density profile and a nonlinear analysis of the orbital 
resonances. Higher resolution simulations clearly illustrate the importance of discrete resonances in 
transporting the angular momentum from the bar to the halo. 

Subject headings: galaxies: structure — galaxies: evolution — galaxies: kinematics and dynamics - 
methods: N-body simulations — cosmology: dark matter 
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1. INTRODUCTION 

More than 2/3 of d isk g ala xies host stel lar bars 
(e.g.. iKnaoen et all 12000: iGrosbol et all 12004 
Ma rinova fc Joged l2007h and evolution of this frac- 
tion with redshift is a matter of an on going debate 
(e.g. , Uogee et all l2004t ISheth et all l2008h . Numerical 
simulations of disk galaxies have shown that bars form 
eithe r as a result of a global gravitational in stability 
(e.g.. lToomrdll98H:ISeilwood fc Wilkinsonlll9p) or thev 
are t r iggered by gal axy interactions (e.g., iByrd et all 
119861: iNoguchil Il987ft and inter a ctions with DM sub- 
struc t ure (e.g., iGauthier et al.l 120061 : iDubinski et al.l 
120081: iRomano-Diaz et al.l l2008f >. A large body of 
theoretical work on the bar instability has examined 
the properties of bars that emerge in initially unstable 
disks in A-body simulations. While these experiments 
explore an idealized picture of bar formation, they 
reveal important aspects of the phenomenology of 
the bar instability, including bar growth within the 
corotation (CR) radius, the vertical buckling instability, 
and the transport of angular momentum through 
gravitational torques from resonant orbits in the outer 
disk and the surrounding dark matter (DM) halo. 
The importance of the resonance nature of angular 
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mo mentum loss by bars an d spira ls was first pointed out 
by iLvnden-Bell fc Kalnaid (|1972h . Angular momentum 
transfer was studied subsequently b oth in ide alized 
models with rigid bar s in li v e halos (Weinberg! 
Hernquist fc Weinberg! Il992t IWeinberg fc Katzj 
a[); and self-consistent A-body simulations 



1985 



2002 



2007aD ; and self-consist ent A-body simulations with 
bar-unstable disks (e.g., ISellwood I Il980l: lAthanassoulal 



1996; 



2003; 



Debattista & Sellwood 1998; Valenzuela & Klypin 



O'Neill fc Dubinskill2003f ). with resonant transfer 



mec hanisms being explored explicitly in some st udies 
(e.g. [A thanassoula 200 21 iHollev-Bockelmann et alJl200r ' 
Martinez- Valpucst a et alj 120061 : iCeverino fc K lypi 



2001 ? These studies have shown that the halo absorbs 
angular momentum from the bar that leads to the 
decline of the bar pattern speed. 

Previous results reveal a close connection be- 
tween numerical bars and observed galactic systems 
in many structural details, including a link be- 
tween the peanut-shaped bulge s and (buckled ) bars 
e.g.. iCombes fc Sanders! [19811 iCombes et al.l Il990t 



iRaha et al.l 119911: iBerentzen et al.1 119981: iPatsis et al 
| 2002t iMartinez-Valpuesta et all 120061 : iDebattista et al 
|2006|) . The ob servational determination of bar pat- 
tern speeds (e.g.,lKentJll987l : lMerrifield fc Kuijkenlll995l : 
ICorsini et al.l [2007) sug gest that stellar bars ar e predom- 
inantly "fast" (but see lRautiainen et al.l (2008) for a dif- 
ferent view) meaning that they are near the maximum 
possible length of the CR radius permitted by the or - 
bital dynamics (jContopoulosI 119801 lAthanassoulal [19921. 
If evolved for too long, the numerical bars can appear 
"slow" with lengths significantly shorter than the CR 
radius and pattern speeds that seem abnormally low 
when compared to observations of real barred galaxies 
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(jDebattista fe Sellwoodlll998ll2000f) . However, at higher 
resolution, even collisionless numerical bars seem to grow 
in length towards their CR radius by capturing disk 
orbit s and so remain "fast" (|Martinez-V alpucst a et all 
2006). Furthermore, the addition of gas may stabilize 
the bar against braking and results in its speedup inste ad 
for prolonged time periods (| Romano- Diaz et al.ll2008[h 

Some studies claim that bars may destroy the cuspy 
profil es of DM halos predicted by the CDM cos molog y 
(e.g., iDubinski fe Carlberd Il99lf iNavarro et ail Il996f ). 
thus alleviating an apparent contradiction between the 
inferred density profiles of DM halos from galactic rota- 
tion curves and this theoretical expectation in some cases 
(I Weinberg fc Katzll2002HHollev-Bockelmann et al.ll2005t 
IWeinberg fc KatT 2007bf K Simulations demonstrating 
cusp destruction use rigid, ellipsoidal bars — their ap- 
plicability to self-consistent dynamical systems is sus- 
pect. Also, there has been some concern about artifi- 
cal m = 1 instabilities arising from usi ng a fixed cen- 
ter in iV-body field expans ion methods (Scllwoodl 120031 : 
iMcMillan fc Dehnenl I2005D . Current studies have ob- 
tained contradictory results on the efficiency of angu- 
lar mom entum transport to the cusp. IWeinberg fc Katzl 
(|2007af) have emphasized the importance of numerical 
resolution and, specifically, of the total particle number 
in simulations. Since the transport of angular momen- 
tum operates mainly through low order resonances be- 
tween the bar pattern speed and halo orbital frequen- 
cies, only a small fraction of the halo mass participates. 
Without adequate particle numbers then, they argue 
that torques associated with resonant populations may 
be under-sampled, leading to a spurious calculation of 
angular momentum tra nsport and, therefo r e, the e volu- 
tion of the bar overall. IWeinberg fc Katzl (|2007bf ) esti- 
mate that at least 10 8 particles and maybe more may 
be necessary to sample the phase-space densely enough 
to converge to the correct answer. Recently, Sellwood 
(2008) has disputed this claim in simulations with rigid 
bars in spherical, isotropic halos with ~ 10 8 particles 
arguing that the resonances are broader than they claim 

In this paper, we address the issue of the numeri- 
cal convergence of bar evolution using a series of N- 
body simulations of the bar instability in a self-consistent 
model galaxy. We analyze bar growth in a bar-unstable 
./V-body disk. In contrast to other work, we em- 
ploy new galactic models based on the methods of 
IWidrow fc Dubinskil (|2005f ). and carry out simulations 
with substantially greater numerical resolution than re- 
ported in the literature. The galaxy is described by a 
well-defined distribution function for an exponential disk 
embedded within a DM halo with an r -1 density cusp, 
based on a truncated Navarro, Frenk & White (1996, 
NFW) profile. These models are formally in dynami- 
cal equilibrium but are bar-unstable. Since they are de- 
fined by a distribution function, their iV-body realiza- 
tions are equivalent, independent of the particle num- 
bers. Hence, this study can probe the effect of numerical 
resolution on collisionless galaxy evolution. Our goal is 
to quantify the behavior of a number of specific param- 
eters describing the bar instability as a function of par- 
ticle number, including the bar strength amplitude, A2, 
as given by the m = 2 Fourier mode, as well as its pat- 
tern speed evolution, angular momentum transport, and 
evolution of the DM density profiles, particularly in the 



region within the halo characteristic NFW scale radius, 
r s . We also perform an orbital spectral analysis of halo 
and disk particles, to quantify the effect of the low order 
resonances responsible for angular momentum transp or t 
(|Athanassoulall2002trMartinez-Valpuesta et al.ll2006| i. 

The plan of the paper is as follows. In §2, we provide 
a description of the galactic models and the TV-body ex- 
periments to study the bar instability. In §3, we present 
results on the bar growth and the evolution of pattern 
speed as a function of numerical resolution. In §4, we 
examine the evolution the DM halo density profile as a 
function of numerical resolution. In §5, we study the low 
order resonances between the bar and the halo particles 
using orbital integrations and spectral analysis and again 
compare results at different resolutions. We also exam- 
ine the details of the evolution of the halo phase space 
density in our highest resolution models. We conclude 
with a discussion of the importance of numerical resolu- 
tion in these experiments and comment on the reliability 
of current work in studies of disk galaxy formation and 
dynamics. 

2. METHODS 

2.1. Initial conditions: An exponential disk with a cuspy 

dark halo 

The main goal of this study is to characterize the bar 
instab ility in terms of mass resol ution. The galaxy mod- 
els of IWidrow fc Dubinskil (|2005f) (WD models herein) 
are ideal for this purpose since they are derived from 
a composite 3-integral distribution function (DF) / = 
fdisk(E,L z ,E z ) + fhaio{E). The disk model has an ex- 
ponential radial profile and sech 2 z vertical profile. The 
disk DF fdjs k is a 3D extension of the 2D function in- 
troduced by[Sh3 (1969) using the vertical energy E z = 
l/2z 2 + &(R, z) — &(R, z = 0) as a n approximate third 
integral (jKuiiken fc Dubinskil fl995h . This DF applies 
in the epicyclic approximation with ur.^.z <C v c and 
so the vertical energy is approximately constant. This 
leads to triaxial velocity ellipsoids in the disk models as 
seen in real spiral galaxies. These models generally pro- 
vide near equilibrium initial conditi ons and show negli- 
gible transient behavior at startup (|Widrow fc Dubinskil 
2005). The halo DF fhalo describes a truncated spheri- 
cal, isotropic NFW model. When the two DFs are com- 
bined, the net halo density profile changes slightly from 
the NFW form and is flattened along the z-axis near 
the center, but preserves the r _1 central cusp. A suit- 
able choice of parameters allows the construction of a 
realistic model of bulgeless spiral galaxy with a cosmo- 
logically inspired DM halo. Since the models are derived 
from a distribution function, particle distributions for 
iV-body experiments can be generated by direct Monte- 
Carlo sampling. 

For the experiments described below, we initially gen- 
erate a model containing 18M disk particles and 100M 
halo particles with both disk and halo particles having 
approximately the same mass. The halo is non-rotating. 
Lower resolution models are generated by subsampling 
this larger model in factors of ten and hence creating a 
sequence of models containing numbers of particles in the 
range 1.18 x 10 4 ~ 8 . One further model is generated with 
a multi-mass DM halo to increase the particle number 
density in the core by another two orders of magnitude. 
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The particle mass is weighted as an approximate step 
function in angular momentum m ~ m(L) such that low 
angular momentum particles near the halo center below a 
characteristic angular momentum L c would have a lower 
mass. The number density at the center of this model is 
more than 100 x greater so the effective particle number 
is » 10 10 for this simulation. We describe the details for 
generating the multi-mass model below. Our highest res- 
olution simulations have large enough particle numbers 
to probe the divergence in n umerical behavior discussed 
by ([Weinberg fe Katzll2007al) . 

Each model is generated and simulated in units with 
G = 1 and physical quantities are of order unity. We 
have designed the model as a proxy for the Milky Way 
without a bulge, so natural units for this comparison are 
L = 10 kpc, M = 10 11 M , V = 207.8 km s" 1 and 
T = 47.2 Myr. By design, the model mass profile closely 
resem bles the one examined bv lMartinez-Valpuesta et al.l 
(2006). Moreover, the central density cusp is better re- 
solved and the initial conditions are in a better equilib- 
rium, since they are sampled from a DF. Throughout this 
paper we present results in physical units. 

The galaxy mass model is presented in Figure[T]as a ro- 
tation curve decomposition. We use an exponential disk 
with radial scale- length 2.85 kpc and an exponential ver- 
tical scalelength of 250 pc and a total mass 5.5 x 10 10 M Q . 
The disk is truncated smoothly at R = 21 kpc equivalent 
to 7.4 scale lengths. The NFW halo scale radius in the 
DF of the WD model is set to r s — 10 kpc but results in 
an effective scale radius of r s = 4.3 kpc as measured by a 
least-squares fit to the density profile. The peak circular 
velocity of the DM halo is v max = 0.77 (160 km s _1 ). 
We note that the smaller scalelength is not due to an 
adiabatic contraction, but is the result of combin ing two 
distribution functions (|Widrow fe Du binski 2005T) — the 
extra concentration of mass from the potential of the disk 
causes the halo potential derived from the NFW DF to 
be more concentrated as well when calculating the self- 
consistent potential for the model. The halo extends to 
a truncation radius of r = 260 kpc and has a total mass 
M = 3.0 units (3.0 x 10 11 M Q ). The final model is a real- 
istic facsimile of an exponential disk galaxy with a cuspy 
DM halo. The square of the radial velocity dispersion 
a 2 R of these models follows the same exponential radial 
decline as the surface density with a\ ~ exp(—R/Rd). 
We choose a central value ctr,o = 104 km s^ 1 , so that 
the Toomre Q is Q = 1.1 at R = 10 kpc. The disk is, 
therefore, relatively cold and responsive. This model is in 
dynamical equilibrium but also is strongly bar-unstable. 
Our analysis focuses on the development of the bar insta- 
bility in simulations of this model with different particle 
numbers. At this point, we also present the final state 
of the mass model after 9.4 Gyrs of dynamical evolution 
for direct comparison to the initial state but defer the 
discussion until later (Fig. [2]). 

2.2. Multi-mass model 

A common way of increasing mass resolution with a 
number of particles is to use a range of masses, assigning 
low mass particles to the center where the action is and 
high mass particles to the periphery ([Sigurdsson et al] 
1995). We therefore build an additional model that 
weights the halo particle mass as a monotonically increas- 
ing function of orbital angular momentum L — r x v|, 
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Fig. 1. — Initial circular velocity curve of the mass model show- 
ing the contributions from the disk and the DM halo. We also plot 
the mean tangential velocity in the disk to show the effect of an 
asymmetric drift on the rotation curve. 
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Fig. 2. — Final circular velocity curve of the barred galaxy 
mass model at t = 9.4 Gyr. We show the contributions from the 
barred disk and the DM halo. The disk component is estimated by 
axisymmetrizing the barred disk and calculating v 3 = RdQ / dR. 



to increase the number density of particles in the region 
where the bar forms and where the low order resonances 
occur. The strategy is to define a mass weighting func- 
tion W(L) such that particles with low angular momen- 
tum and orbits with small pericenters also have small 
mass, while those with large angular momentum and 
pericenters beyond the edge of the disk have a higher 
mass. The halo DF is normalized by this weighting func- 
tion, so that the number density of particles derived from 
Monte Carlo sampling will be larger for smaller values of 
L. In this way, the probability of selecting a particle 
with smaller L is greater than with larger L. The biased 
number density is then corrected to represent the model 
with the original DF by multiplying the particle mass by 
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W(L). The weighting function is normalized so that the 
mass of a particle in the initial distribution is given by 

M halo W(Li) 

The choice of the functional form of W(L) is arbitrary 
at some level according to the needs of the problem but 
in our case it should be monotonically increasing with L. 
We use the step-like weighting function in L 

where L c is a characteristic angular momentum for the 
step, a > is an exponent and Wi is the asymptotic 
value of weighting function for large angular momentum. 
When W is plotted versus log L it takes the form of a step 
function where the steepness of the transition at log L c 
depends on the choice of a. In practice, we truncate the 
function at minimum and maximum values of L at L m i„ 
and L max and set the weight to the value at these limits 
beyond the endpoints. 

After some experimentation, our final choices for these 
parameters are Wi = 10 4 , L min = 1CT 3 , L c = 3, L max = 
7, and a = 0.9. The choice of L c corresponds to particles 
moving at the circular velocity at a radius of R — 4.1 (41 
kpc) about twice the radius of the disk. The choices of 
L rn i n and L max limit the dynamic range of masses to 
about 600 with the least massive particles weighing in at 
0.5% the equivalent mass for a single-mass model and the 
most massive particle weighing in at 3 x the equivalent 
mass. For comparison, the single-mass particle in the 
N = 10 8 halo is 3 x 10 3 M Q , while in the multi-mass 
model, the particle masses range from 16 M for small 
L to 10 4 M for the most massive particles in outskirts 
of the halo. 

We plot the ratio of the particle number density in 
the multi-mass model to the equal mass particle number 
density in Figure [3] The number density is about 200 x 
greater within the central 100 pc of the model and about 
10 x at R = 1 kpc. 

2.3. Simulations 

We simulate t hese models using a parallelized treecode 
iDubinskil (|1996l) for 200 time units (9.4 Gyr), permitting 
us to see the development of the bar instability through 
various phases roughly over a Hubble time. We soften 
gravity with a Plummer model kernel and vary the soft- 
ening length e according to the particle numbers of the 
simulation roughly in proportion to N^ 1 / 3 . The median 
force errors are 0.1% for the chosen treecode parameters. 
Simulation parameters are given in Table 1. 

A constant timestep is used for all of the simulations 
(see Table 1). The circular orbital period in the mass 
model at the smallest softening radius of e = 10 pc 
is about 15 Myr and so is resolved by 30 timesteps. 
Plummer softening smooths gravity over a few softening 
lengths so the smallest resolved radius of these simula- 
tions is « 3e. We find that total binding energy is typ- 
ically conserved to within 0.2% and angular momentum 
is conserved to within 1% over the course of the runs. 
Each simulation, with the exception of the multi-mass 
case, is repeated twice with a different random realiza- 
tion to explore statistical variance in the growth of the 
bar mode. 



r (kpc) 

Fig. 3. — The ratio of the number density of the multi-mass 
100M particle halo to the single-mass 100M halo. The distribu- 
tion function is sampled such that particle mass is weighted by 
a smoothed step function of total angular momentum. Particles 
with low angular momentum have small mass and those with high 
angular momentum have low mass (see text). The particle density 
is more than 100 times higher within 0.1 kpc and at least 10 times 
higher within 1 kpc. For R > 10 kpc the number density drops 
gradually to about half the single-mass case. The effective numeri- 
cal resolution at the center of simulation is therefore Nh ~ 10 9 ~ 10 . 



TABLE 1 

Simulation parameters 



Model 


N h 




N d 


e (pc) 


St (kyr) 




mlOK 


10 4 


1 


8 x 10 3 


200 


470 


20000 


mlOOK 


10 5 


1 


8 x 10 4 


100 


470 


20000 


mlM 


10 (i 


1 


8 x 10 5 


50 


470 


20000 


mlOM 


10 7 


1 


8 x 10 6 


20 


470 


20000 


mlOOM 


10* 


1 


8 x 10 7 


10 


470 


20000 


mmlOOM 


10* 


1 


8 x 10 7 


10 


235 


40000 



Note. 
model. 



The model mmlOOM 



the multi-mass 



Figure 2] shows an animation 5 of the evolution of 
the disks in six models in face-on and edge-on views. 
The lowest resolution model mlOK demonstrates how 
insufficient particle numbers can lead to spurious re- 
sults. A bar develops immediately in the 1.8K particle 
disk but devolves into a compact rapidly tumbling ob- 
ject. In retrospect, early galaxy formation simulations 
that introduced the angular momentum problem (e.g., 
iNavarro &: Stei nmetz 2000) only contained 2K particles, 
so part of the problem may have arisen from exceed- 
ingly noisy evolution of a bar mode. The mlOOK model 
with an 18K particle disk still appears noisy, though the 
buckling instability is clearly visible. The disk is visibly 
thicker than the higher resolution models, however, and 
the bar is not as pronounced. Disk heating by bombard- 
ment of halo particles is a problem. The time of onset 
of the bar instability is delayed as N increases, reflect- 
ing the effect of Poisson noise. Since the bar instability 
grows exponentially from density fluctuations in the ini- 

Quicktime animations are available at the website 
www.cita.utoronto.ca/~dubinski/BarsInCuspyHalos/ 
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tial conditions, larger N simulations will have smaller 
initial amplitudes and, therefore, longer times to satu- 
rate. 

Figure [5] displays the evolution of the multi-mass halo 
model with the 18M particle disk close-up and two per- 
pendicular edge-on views simultaneously. Figure [6] refers 
to the face-on view in a frame co-rotating with the bar 
to emphasize the growth of the bar mode. This model 
starts very quietly and there is little visible structure un- 
til t Ri 1 Gyr when the bar begins to emerge. The bar 
grows from the inside out, gradually increasing in length 
until it reaches a maximum length at nearly the CR ra- 
dius around t = 2 Gyr. At this time, the bar also ex- 
cites a prominent, bi-symmetric spiral structure. After 
saturation, the bar re-structures itself, becoming more 
centrally concentrated and weakening, as it settles into 
a quasi-steady state. After settling, the pattern speed 
begins to decline and the bar's length increases slowly, 
since the CR radius is increasing and the bar can capture 
additional orbits in the disk. The other notable event 
is the vertical buckling instability that occurs around 
t = 3.5 Gyr creating a characteristic X-shaped structure 
as various families of orbits establish themselves causing 
the bar to thicken vertically. By the end of the simula- 
tion, the inner bar transforms into a peanut-shaped bulge 
though it still is obviously elliptical in the face-on view. 

Figure [21 shows the rotation curve decomposition for 
the model at the final time t = 9.4 Gyr. We rotate parti- 
cles in the final barred disks to random angles <f> to make 
the disk potential axisymmetic and then compute the 
disk rotation curve component through v 2 d = Rd^/dR. 
The halo rotation curve component is estimated from the 
spherically-averaged density profile of the dark matter. 
The collapse of the bar leads to a concentrated bulge- 
like component and the rotation curve flattens slightly 
creating a galaxy model that more closely resembles real 
systems. The halo profile at large radii does not change 
a lot but we will see in further analysis discussed below 
that there is slight increase in the central density. In 
this barred galaxy model, both the stars and dark mat- 
ter have comparable contributions to the rotation curve 
in the inner regions. 

In the next section, we quantify these various effects 
and look for differences in resolution with the hope of 
finding numerical convergence in physical behavior. 

3. BAR GROWTH AND PATTERN SPEED EVOLUTION 

The growth of the bar instability is measured by the 
bar strength with the m = 2 Fourier amplitude of the 
surface density, A 2 , given by: 

1 N 

A 2 = ~Y^ m i exp(2i(f>j); R < R c (3) 
i=i 

where the summation is performed over a list of particles 
with masses m, at angle 4>i in the x — y plane within some 
cut-off radius R c . The normalized amplitude \A%\ versus 
time measures the growth rate of the bar instability and 
the phase angle 4> = 0.5 tan _1 [Im(A2)/Re(A2)] with time 
permits measurement of the pattern speed by numerical 
differentiation. 

All models were simulated for 200 time units (9.4 Gyr) 
and every 10 (or 20 in the multi-mass case) timesteps 
a face-on surface density image was generated from the 



particle distribution resulting in a sequence of 2000 im- 
ages. These image sequences were analyzed to deter- 
mine A2 by summing over pixels rather than particles in 
equation [3j The amplitudes and phase angles were then 
tabulated as a function of time to determine the rate 
of growth of the bar and the pattern speed evolution. 
Pattern speed is estimated by simply differencing angles 
in subsequent pairs of phase angles and dividing by the 
time interval. In practice, we only use values every 50th 
snapshot to smooth out the noise in these parameters 
introduced by limited numbers of particles. We show be- 
low that higher resolution simulations produce smoother 
curves of bar growth and pattern speed evolution. 

3.1. Numerical Accuracy 

We first discuss the behavior of t he bar instab i lity as 
a function of integration timestep. iKlypin et al.l (|2008[ ) 
have claimed that very small timesteps are necessary to 
resolve the dynamics of bars because of the possible de- 
velopment of cuspy density profiles in the forming bulge- 
bar system. Our simulations use a single timestep chosen 
to resolve the smallest dynamical timescale in the model. 
Multiple timestepping schemes often use the criterion, 
St — (2.8e/|g|) 1 / 2 ?7 where e is the Plummer softening or 
equivalent, g is the acceleration and 77 is a free parame- 
ter us ually chosen with a recommended value of 77 = 0.2 
(e.g.. ISpringel et aT1l2001h . For density laws following 
p ~ the central acceleration is constant. The high- 
est value of the acceleration in our galaxy model occurs 
in the center, and using it with rj = 0.2 we arrive at 
St = 0.01 (470 kyr) for e = 50 pc and St = 0.004 (190 
kyr) for e = 10 pc. Plummer softening of course reduces 
the maximum value near the center and it formally falls 
to zero at r = 0. We see below that there is only a 
modest increase in the central density evolution, so the 
maximum value of g does not change by much over the 
course of the run. 

The smallest orbital period is ~ 15 Myr for an or- 
bit with R ~ e and our chosen timestep is St = 470 
kyr, so these orbits are resolved with approximately 30 
timesteps. For our highest resolution simulation, we use 
St = 235 kyr to account for the smaller softening radius. 
The fraction of particles with orbital periods less than 
20 Myr is approximately 0.1% based on a analysis of the 
radial frequency of 100K testparticle orbits sampled from 
the halo integrated within the rigid potential of our mass 
model. If the timestep is too large, orbits near the center 
will be unstable and create an artificial constant density 
core. Another possible problem occurs for highly radial 
orbits with longer periods that pass close to the central 
cusp. Our orbital analysis showed that approximately 
0.13% of orbits change binding energy by more than 1% 
over a 9.4 Gyr integration with St — 470 kyr. All of these 
orbits had small pericentric radii ~ 100 pc. We therefore 
expect a small fraction of highly radial orbits to diffuse 
artificially through energy space. We demonstrate here 
that the single timesteps of St = 235 and 470 kyr are 
sufficiently small to resolve the dynamics for our choices 
of the Plummer softening radius. 

To test for numerical convergence, we have re-run the 
model mlM using single timesteps over the range of St = 
15 - 940 kyr for a time of 4.7 Gyr. The model with 
St = 15 kyr required 320K single steps. This galaxy 
model has Nd — 180K and Nh — 1M and a Plummer 
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Fig. 4. — A comparison of the evolution of the bar instability in 6 simulations with increasing particle number N. The formation of a 
bar is delayed for simulations with larger N since the Poisson seed noise has a lower amplitude and it takes longer for the instability to 
grow in the linear regime. The lowest resolution simulations suffer from heating while the general behavior converges at higher resolution 
for TV > 10 6 (see Video 1) 




Fig. 5. — Evolution of the multi-mass model in the inertial frame showing the face-on view and two perpendicular edge-on views. The 
bar grows from the inside out first evolving into a thin bar extending to the co-rotation radius and then settling down into a less elongated 
ellipsoid. The buckling instability vertically thickens the bar into a peanut-shaped bulge at later times. The bar grows in length as angular 
momentum is lost to the halo and new orbits are captured with the co-rotation radius (see Video 2). 



7 





0.002 0.004 0.006 0.008 0.01 

(5a/ a 



Fig. 7. — Relative acceleration errors for the parallel treecode for 
runs with N = 1.18M particles. Errors are estimated by comparing 
acccelerations computed with our preferred treecode opening angle 
parameter 9 = 0.9 with quadrupole order corrections to the exact 
values determined from a direct calculation. The mean relative er- 
ror is 0.13% with a median value of 0.085%. Note the opening angle 
criterion for the parallel treecode is more conservative than stan- 
dard definitions and so a larger value of 6 still results in relatively 
small acceleration errors (Dubinski 1996). 

softening length of e = 50 pc. We examine different 
metrics of the system evolution including bar growth, 
pattern speed evolution and the final density profile of 
both the stars and DM all as a function of timestep. 

3.1.1. Acceleration Errors 

We first comment on the accuracy of the accelera- 
tions determined using the parallel treecode (Dubinski 
1996). Normally force accuracy is not discussed despite 
a large variety of algorithms used to compute gravita- 
tional forces. We present our errors here so that other 
researchers may compare to their own standards of nu- 
merical accuracy. Figure [7] shows the distribution of rel- 
ative acceleration errors for our preferred treecode pa- 
rameters. We use an opening angle tolerance 6 = 0.9 
with quadrupole corrections using a more conservative 
cell opening criterion than normally described that gives 
more accurate acceleration values for a given 9 than stan- 
dard definitions (Dubinski 1996). Errors are determined 
by comparing accelerations from the treecode method to 
a direct force calculation. The median and mean rel- 
ative acceleration errors are 0.085% and 0.13% respec- 
tively with 99.7% (3er limit) of acceleration errors less 
than 0.7%. 

3.1.2. Energy Conservation 

The evolution of the error in total binding energy of an 
N-body system is a useful indicator of the fidelity of the 
results and can reveal potential problems with the inte- 
gration scheme or choice of timestep. Figure |5] shows the 
change in total binding energy as function of timestep. 
The largest timestep of St — 940 kyr shows a strong sys- 
tematic drift in energy reflecting the inadequate timestep 
resolution for a significant fraction of orbits. There is a 
smaller drift in the energy with a relatively small error 
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Fig. 8. — Total energy errors for runs with different time-steps 
St. The simulation with St = 940 kyr shows a systematic drift 
due to inadequate numbers of timesteps to follow orbits within the 
core. There is a lesser drift for the timestep St = 470 kyr but the 
error has only grown to 0.1% by the end of the run. All timesteps 
with St < 470 kyr show very little drift. 

of 0.1% over 4.7 Gyr with our main choice of St = 470 
kyr and clear convergence with no systematic effects with 
St < 235 kyr. We, therefore, conclude that St = 470 kyr 
is adequate for our models. We show below that there 
are no substantial differences between various metrics of 
the properties of the bar and halo when using timesteps 
with S < 470 kyr. 

3.1.3. Bar Evolution versus Timestep 

We measured both bar growth and pattern speed evo- 
lution as a function of timestep in the mlM model with 
N d = 180K and N h = 1M. Figure M presents the evolu- 
tion of the bar growth parameter | ^4.2 1 measured within 
R < 5 kpc versus timestep. During the linear growth 
phase of the bar instability, all simulations track one 
another very closely. However, after the bar instabil- 
ity saturates around t ~ 1 Gyr the behavior is quite 
variable and erratic for different choices of the timestep. 
The time of bar buckling shown by the sudden secondary 
drop in \A2\ changes with different timesteps and lies in 
the range t — 1.8 — 2.5 Gyr. There is no monotonic 
trend with timestep. The range of variability is the same 
as our study of independent random realizations in § 13.21 
The root cause of this behavior is probably the dynamical 
chaos inherent to this late evolution of the bar instability. 
The detailed TV-body solutions for individual particles 
diverge exponentially for different choices of integration 
step in the nonlinear regime of dynamical evolution. De- 
spite this divergence on the individual particle level, the 
global properties of the resulting bar are similar as we 
shall see. 

An analysis of the pattern speed evolution shows con- 
sistent results for all timesteps (Fig. [10]). The agree- 
ment in the linear regime evolution until the bar in- 
stability saturates at t ~ 1 Gyr is very close, after 
which the detailed evolution show differences. There is 
a 2 — 4 km s _1 kpc -1 scatter in the pattern speed at 
any given time but the general declining trend is the 



9 




012345 01234 
t (Gyr) t (Gyr) 



Fig. 9. — Evolution of the Fourier component, A 2 for stars with 
R < 5 kpc for a single model with N d = 180K and N h = 1M 
particles with global timesteps spanning the range of St = 15 kyr 
to 940 kyr. The linear growth phase of the bar is almost identical 
until the bar instability saturates at t ~ 1 Gyr. The subsequent 
nonlinear evolution shows a wide range of behavior for different 
timesteps with no monotonic trend. The nonlinear phase of the 
bar instability involves chaotic orbits and so the slight variations 
introduced by the round-off error of different discrete timesteps 
lead to divergent evolutionary behavior. The main manifestation 
of this chaos are different times for the onset of the buckling in- 
stability ranging from 1.8-2.5 Gyr having no dependence on the 
chosen timestep. Nevertheless, the behavior is qualitatively simi- 
lar after the buckling instability with a steady rise of \A%\ at late 
times as the bar lengthens. 

same over the course of the run. The observed scatter 
is consistent with the same scatter seen in different ran- 
dom realizations (Fig. I13[) . The mean and variance of 
the pattern speed at t — 4.7 Gyr for all timestep runs is 
£!b = 16.6 ± 0.7 km s _1 kpc -1 . We conclude that our 
choice of timestep leads to a consistent evolution of the 
bar pattern speed. 

3.1.4. Stellar and Halo Central Density versus Timestep 

As a final metric of the accuracy of the simulations 
versus timestep, we measured the spherically averaged 
density profile of the stars and DM at the last snapshot 
at t = 4.7 Gyr. At this time, the bar has buckled and 
has formed a concentrated bulge-like object within the 
halo that has become more dense itself in response to 
this new bulge (see below) . Figure [TT] shows the stel- 
lar density profile within r < 1 kpc of the center for 
runs with different timesteps. The density profiles for 
the bar/bulge are consistent within the error bars for 
St < 470 kyr. There is some random scatter in the inner 
radial bins since there are only a few hundred particles at 
these small radii. The model with St — 940 kyr forms a 
core with constant density within r < 200 pc though the 
density is only 0.3 dex (about 2x) smaller than the den- 
sity in the first radial bin of the smaller timestep runs. 
The runs with St < 470 kyr agree within ±0.1 dex for 
r < 200 pc and much of that error is due to small particle 
numbers with ~ 10 2 particles per bin. 

We also measured the spherically-averaged density pro- 
file of the DM at t = 4.7 Gyr (FigfT2j). We have approx- 
imately 6 x as many particles per bin and so the random 



Fig. 10. — Evolution of the pattern speed Qt, for models with 
JVd = 180K and = 1M and timesteps spanning the range of 
St = 15 kyr to 940 kyr. During the linear growth phase of the 
bar until t s=s 1 Gyr, the pattern speed evolution is almost the 
same. Once the bar becomes nonlinear, there is a small scatter 
in the detailed behavior of the pattern speed with a variation of 
2 — 4 km s — 1 kpc -1 at any given time. By the end of the runs, 
the results converge with the mean and variance of the pattern 
speed £7b = 16.6 ± 0.7 km s _1 kpc -1 at t = 4.7 Gyr. There is no 
strong dependence of the pattern speed evolution on the choice of 
timestep. 

errors are smaller. The DM density profiles are consistent 
for St < 470 kyr suggesting we have adequate time reso- 
lution for the halo density evolution. Again, we see the 
development of an artificial constant density core in the 
simulation with St = 940 kyr. This timestep is clearly too 
large and does not adequately follow short period orbits 
in the core. However, simulations with timesteps smaller 
than St < 470 kyr adequately follow the dynamics of the 
evolution of the DM halo. 

In summary, we have presented the force accuracy 
and total energy evolution of our simulations with dif- 
ferent timesteps. We have also shown that our results 
converge experimentally for St < 470 kyr according to 
different metrics of the bar evolution including pattern 
speed evolution and stellar and dark matter density cen- 
tral density profile. We note that there is a random be- 
havior for the time of onset of the buckling instability 
for different choices of timestep which probably results 
from the chaotic nature of this dynamical system. (This 
was shown explicitly by Martinez- Valpuesta & Shlosman 
2004.) If the timestep is too large, the main effect is 
to create an artificial constant density core. Particles 
with short orbital periods are numerically unstable and 
are scattered out of the center creating the core. In the 
subsequent analysis, we show that the central density 
continues to increase at smaller radii with higher mass 
resolution. If our timestep was too large, one might ex- 
pect instead to see the onset of a artificial constant den- 
sity core of a fixed radius set by the timestep and inde- 
pendent of the mass resolution. We do not observe this 
behavior. We also do not see a sudden change in behav- 
ior of the pattern speed evolu tion at a critical timestep 
as seen bv lKlypin et al.l (|2008f ). We, therefore, conclude 
that we have adequate time resolution to follow the dy- 
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Fig. 11. — The spherically averaged density profile of the stellar 
component that includes the buckled bar and disk at t = 4.7 Gyr 
for runs with different time-steps. The error bars are 1— a estimates 
of the \fN Poisson error in the density due to discrete sampling 
e.g., the inner most bins contain ~ 100 particles so the 1 — a error 
in density is about 10%. For time-steps with St < 470 kyr, the 
density profiles are consistent within the random errors. The run 
with St = 940 kyr shows the formation of an artificial core due to 
an insufficient number of time steps to follow orbits within r ~ 100 
pc. 
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Fig. 13. — Evolution of the Fourier component, A% for stars 
with R < 0.5 for 10 models with N d = 180K and N h = 1M 
generated with different initial random seeds. There is a large 
variation in evolution of Ai during the formation of the bar over 
the time interval t = 20 — 70 reflecting detailed differences in the 
Poisson noise in different random realizations. The plot reveals the 
approximate scatter in evolutionary behavior expected for different 
runs. 
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Fig. 12. — The spherically averaged density profile of the dark 
matter halo at t = 4.7 Gyr for runs with different time steps. The 
error bars are 1 — a estimates of the s/N Poisson error in the density 
due to discrete sampling. For time-steps with St < 470 kyr, the 
density profiles are consistent within the random errors. The run 
with St = 940 kyr shows the formation of an artificial core due to 
an insufficient number of time steps to follow orbits within r ~ 100 
pc. 

namical evolution of this system all the way down to the 
radius where Plummer softening dominates. 

3.2. Models at fixed resolution 

Before presenting results on the bar and pattern speed 
evolution versus mass resolution, it is instructive to un- 



derstand the variance expected for runs at a fixed reso- 
lution. The seed of both spiral and bar instabilities in 
./V-body simulations is the Poisson noise in the discrete 
particle distribution of the disk and halo. We, there- 
fore, expect some variation in the detailed behavior of 
the growth of the bar instability in different random re- 
alizations and we quantify it here. 

We build ten galaxy models with 1M halos particles 
and 180K disk particles independently from different 
Monte-Carlo samplings of the galaxy model DF by using 
a different initial seed for the random number generator. 
We measure |^4.2 1 within a radius R < 0.5 units (5 kpc) 
which is within the eventual co-rotation radius of the bar. 
Figure [T3l shows the evolution of the bar strength for the 
10 runs. The detailed behavior varies significantly for the 
different runs with the minimum and maximum values of 
| v4_2 1 that varying by ±0.05 during the bar growth phase 
between t = 1 — 3.3 Gyr and final values ranging from 
0.35-0.40 at t = 9.4 Gyr. While the runs differ in detail 
there is still a generic behavior with the bar growth sat- 
urating around at IA2I ~ 0.5 and then going through an 
oscillation before settling down to a value near \A 2 \ ~ 0.3 
around t = 3.3 Gyr. The bar then grows slowly increas- 
ing in length as the pattern speed declines. 

Fig. [14] shows the pattern speed evolution for the 
same 10 runs at fixed resolution. The behavior is similar 
with the bar starting out with fib s=s 35 km s -1 kpc -1 
declining to a value between 12 — 14 km s _1 kpc -1 . The 
pattern speed evolution is consistent at the 10% level 
despite the different histories of the bar growth as quan- 
tified by \A 2 \. 

We, therefore, expect the minimum and maximum 
values of | ^4.2 1 to vary by about 0.05 units between 
models and pattern speeds to vary conservatively by 
±2 km s _1 kpc -1 for stochastic reasons alone. 
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Fig. 14. — Evolution of the pattern speed, f2f,, for 10 models with 
= 180iC and = 1M generated with different initial random 
seeds. The pattern speed is measured by creating a time series of 
the phase angle of the A2 component for stars with R < 0.5. The 
pattern speed decays after the bar forms as angular momentum 
is transferred to the dark halo through dynamical friction. While 
the decay rate is similar, again there is scatter due to statistical 
variation of the Poisson noise in the initial conditions. 

3.3. Models with increasing mass resolution 

After quantifying the effects of temporal resolution and 
stochasticity in models at fixed resolution, we go on to 
examine models of increasing mass resolution with halos 
containing from 10 5 to 10 s particles and the multi-mass 
model with an effective resolution of 10 10 particles. Our 
goal here is to measure carefully bar growth and pattern 
speed and to single out any differences that are incon- 
sistent with the expected statistical variance. We have 
simulated 2 models at each resolution in Table 1 with the 
exception of the multi-mass case where we did only one 
model. 

Figure [TBI shows the bar growth for all resolutions plot- 
ted as In I A2 versus time to emphasize the growth of the 
instability through the linear regime. Spiral and bar in- 
stabilities grow from seed density fluctuations in Pois- 
son noise thro ugh the swing amplification mechanism 
(|Toomrdfl98l . In the linear regime, the fluctuations 
grow exponentially and so ln|A2| is roughly linear in 
time. The dashed line to the right is parallel to the model 
growth rates and corresponds to exponential growth with 
a timescale of r = 370 Myr. 

Once the perturbation goes non-linear, | ^4_2 1 reaches 
a maximum value and then oscillates until reaching a 
steady state as the bar settles into a quasi-equilibrium. 
All models show a gradual linear rise of In after 
reaching equilibrium but are noticeably offset in the sat- 
uration time when going to higher resolution. Since the 
seed perturbations arise from Poisson noise, the ampli- 
tude of perturbations <5 varies as A -1 / 2 , so the ratio 
of amplitudes in two different simulations is Si/ So = 
(A1/A0) -1 / 2 . In the linear regime, S ~ exp(t/r), so 
the time delay between growing perturbations to reach 
the same amplitude is St « r ^(Ai/Ao) 1 / 2 . Simula- 
tions with a factor of 10 more particles will, therefore, 
be delayed in saturating by a time interval given by 
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Fig. 15. — Initial growth of the bar strength IA2I for 
stars with R < 0.5 for two model sequences using = 
18K, 180K, 1.8M, 18M with = WOK, IM, 10Af, 100M respec- 
tively. The In I A2 1 grows approximately linearly with time indepen- 
dent of the choice of and showing the exponential growth 
of the bar mode. The dashed line shows an exponential timescale 
that is approximately r = 370 Myr. Since the bar grows from the 
Poisson noise within the disk then we expect the noise amplitude 
to be proportional to A r_1 / 2 . Based on exponential growth of the 
bar mode, we expect the time to saturation of \Aj\ to be delayed by 
roughly St ~ Th^iVi/iVrj) 1 / 2 , e.g., a factor of 10 change in particle 
numbers leads to a delay 5t s=s 9. The difference in saturation times 
of \A2 1 between the various simulations are roughly consistent with 
this estimate though there is some variation. 

St ~ rlnl0 1/2 w l.lr. With r w 370 Myr, we expect 
a time delay of approximately 400 Myr between simula- 
tions differing by a factor of 10 in numbers of particles. If 
we select the time when \A 2 \ reaches a maximum as a ref- 
erence time when the bar saturates and the linear regime 
ends, we can estimate the time delay between simulations 
directly. Using the 1M particle run as a zero point, we 
find delay times of St — 280 Myr for 10M particle mod- 
els and St = 600 - 700 Myr for the 100M models and 
St = 950 Myr for the multi-mass 100M model. The noise 
characteristics of the multi-mass model are more com- 
plicated than the simple ideas discussed here and vary 
across the model but the onset of the bar instability is 
nonetheless delayed further because of quieter initial con- 
ditions. These values are slightly smaller than expected 
but are in reasonable agreement with the estimated de- 
lays from considerations of the growth of Poisson fluc- 
tuations. This analysis emphasizes that the spiral and 
bar instabilities that arise in N-body simulations of disks 
are wholly dependent on the initial Poisson noise. In 
the future, with simulations using more than 10M disk 
particles it makes sense to control the properties of the 
noise both in amplitude and power spectrum as done in 
cosmological simulations. 

Figure [in] and [T7] show the evolution of the bar strength 
when we account for the time delays and allow a compar- 
ison of the early and late time evolution. When synchro- 
nized this way, the linear growth phase is readily appar- 
ent in the evolution of In |^4.2 1 in Fig. [TH] In Fig. [T71 the 
plot of the evolution of \A 2 | reveals the details of the non- 
linear evolution of the bar. The bar strength saturates 
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Fig. 16. — Same as Fig. 1 15 1 except the curves have been shifted 
in time so that the linear regime growth phases overlap with the 
mlM model according to the measured time delays 



at a maximum value followed by an oscillation through 
a minimum and then a slow rise to the end of the simu- 
lation. For the most part, the range of behavior between 
different resolutions is consistent with our expectations of 
variance from our study of 10 simulations at fixed resolu- 
tion. However, the highest-resolution multi-mass model 
dips to a lowest minimum value of \A%\ and takes more 
time to grow in the later phase. This difference does 
lie within the range of stochastic behavior but still ap- 
pears slightly anomalous. We will find below that the 
rate of angular momentum transfer between the bar and 
the halo is slightly slower for the multi-mass run. The 
multi-mass run seems to transport about 10% less angu- 
lar momentum from the bar to the halo than the other 
runs and this could account for the different behavior. 

Finally, we compare the pattern speed evolution of 
simulations at different resolutions. Figure [TS] shows 
the pattern speed versus time for all simulations where 
again for a proper comparison we have synchronized 
the various runs to the time of the first peak in | 1 
as before. The decline of the pattern speed is simi- 
lar for all resolutions with bars initially forming with 
ilb ~ 35 km s _1 kpc -1 and ending with a value around 
fib w 12 — 14 km s _1 kpc -1 . The range of curves is 
again consistent with the scatter seen in the fixed res- 
olution study. The highest resolution runs with 100M 
halo particles in both the single mass and multi-mass 
case show an apparent oscillation in during the de- 
cline. The frequency of this oscillation is approximately 
half of the pattern speed Slf, itself. The source of the 
oscillation is not obvious. We initially speculated that 
interference from spiral patterns beyond the end of the 
bar rotating at a different pattern speed may have altered 
the measurement of A 2 within R < 0.5. However, when 
the pattern speed is derived from A 2 measured within 
R < 0.25 out of influence of spirals the oscillations per- 
sist at the same frequency. These oscillations may result 
from uneven bar growth (i.e., variations in length) by 
trapping of the disk or bits by the bar or from a non- 
linear mode coupling (jMartinez-Valpuesta et ail 120061 : 
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Fig. 17. — Evolution of the bar strength \A$\ for stars with R < 
0.5 for two model sequences using N d = 18K, 180 K, 1.8M, 18M 
with N h = 100K, 1M, 10M, 100M respectively and the multi-mass 
model plotted versus linear time. The curves have been synchro- 
nized to the time of maximum bar extent. This plot emphasizes 
the variance in behavior after the bar instability goes nonlinear. 
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Fig. 18. — Evolution of the pattern speed H;, for two 
model sequences using N d = 18K,180K,1.8M,18M with N h = 
100K, 1M, 10M, 100M respectively (R < 0.5). Also, shown is one 
model with N d = 18M and jV/, = 100M with a multi-mass halo 
that increases the particle number density near the center of the 
disk. The curves have been shifted in time so that the bar growth 
evolution is coincident with the mlM model. The decline in pattern 
speed at different resolutions is similar though there the multi-mass 
model does not decay as quickly and has a slightly larger pattern 
speed at the last simulated point. The 100M particle simulations 
also show a modulation of the pattern speed that indicates more 
subtle dynamical effects revealed by higher resolution. 

iMartinez- Valpuestall20M ) . 

In summary, the bar develops from Poisson noise in 
the disks in a similar way for simulations with Nh > 10 6 . 
The time delay in the growth to the nonlinear phase for 
larger Nh are the result of smaller amplitude Poisson 
fluctuations that seed the bar at higher resolution. The 
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variation in the behavior of the different runs is consis- 
tent with the variance introduced from different random 
realizations of the models. The bar pattern speed decays 
at a similar rate over the course of the run for resolu- 
tions again with N > 10 6 though the higher resolution 
runs decay to a final value that is approximately 10% 
larger. There is no dramatic change in dynamical evo- 
lution of gross physical properties of the bar as we ap- 
proach Nh = 10 8 suggesting the models are converging 
to the correct physical behavior. 

3.4. A Fast Bar 

Orbital dynamics permits a bar of length q& to ex - 
tend as far as the CR radius Dl (jContopoulosI Il980h . 
But the developing chaos between the Ultra-Harmonic 
resonance (UHR) and the CR limits the bar length 
to within the UHR, especially in stronger bars. The 
dimensionless ratio 1Z = D^/at, is an indicator of a 
bar's dynamical state and galaxies with observed or 
inferred pattern speeds have 1Z = 1.2 ± 0.2 (e.g., 
lAthanassouh]|1992t[De"battista & Sellwoodlll998t ). Bars 
emerging from the disk instability are usually born 
with 7Z « 1 and this value gradually increases as the 
bar settles into equilibrium and loses angular momen- 
tum to the halo through dynamical friction. Dur- 
ing buckling th e bar shortens dramatically for some 
period of time (IMartinez-Valpuesta fc Shlosmanl 120041 : 



IMartinez-Valpuesta et alj I2006D and afterwards gradu- 
ally lengthens. However, the CR radius also increases 
in response to the change in potential of both the outer 
disk and DM halo as they absorb angular momentum 
from the bar and respond to the changing mass profile 
of the disk. iDebattista &; Sellwoodl (|2000l ) have shown 
that in many models with dense halos, bars are slowed 
down considerably and end up with values of 1Z > 2. 
Bars are, designated as "fast" if 1 < 7Z < 1.4 or "slow" 
for 1Z > 1.4 with all barred galaxies with determined or 
inferred pattern speeds being "fast" by this definition. 

The bar that forms in the model described here is 
"fast" with 1Z < 1.4 after reaching a quasi-equilibrium 
after buckling. The r e sult is in agreement with 
IMartinez-Valpuesta et all (|2006f ) who used a similar 
galactic model and determined the bar size by means 
of the last stable orbit supporting it. At various times 
after the bar forms, we have determined the curve of the 
circular orbital frequency f2(i?) by computing the aver- 
age of the radial acceleration d$ / dR on points on rings of 
different radii in the midplane of the disk. In this way, we 
average out the asymmetry in the potential introduced 
by the bar. The CR radius is then found by reading off 
the radius corresponding to f2(Dr,) = Qb at the given 
time. To determine the bar length we fit elliptical 
contours to the surface density profile and look for a sud- 
den transition in the value of the axis ratio q = b/a and 
the position angle of the isodensity contours. In most 
cases, the transition is sudden, jumping from q = 0.4 to 
q = 0.9 over a radial interval of 1 kpc. We therefore can 
determine with an accuracy of ±0.5 kpc. Figure [19] 
shows the isodensity contours overlayed with corotation 
radius and elliptical contour for the chosen bar radius at 
t = 9.4 Gyr the final time in the simulation for the high- 
est resolution model mmlOOM. Even at this time, the 
bar nearly extends to the CR radius and 1Z = 1.2 ± 0.05. 
Figure [20l shows the evolution of 1Z from t — 4.7 — 9.4 Gyr 
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Fig. 19. — Surface brightness contours of the multi-mass model 
at the last snapshot at t = 9.4 Gyr overlayed with the best fit 
ellipse to the central bar and co-rotation radius. Even at this late 
time in the evolution, the bar extends to the co-rotation radius. 
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Fig. 20. — Evolution of the CR radius to bar length ratio H 
for the multi-mass model over the last half of the simulation. The 
value of H hovers around 1.1 indicating a fast bar. 



starting with the time when it just has settled into equi- 
librium until the end of the run. During this time, the 
bar ratio 1Z maintains a value less than 1.4 and so would 
be classified as a "fast" bar and remain consistent with 
the observed barred galaxies. Despite a 3-fold drop in 
the pattern speed the bar length and the galaxy poten- 
tial r eadjusts to keep 1Z near 1. (|Debattista fc Sellwoodl 
12000ft found that models with V% iak /V% alo = 1 do end up 
with fast bars and indeed our model is consistent with 
that value (Fig. [p. 

Finally, we estimate the bar mass and shape for com- 
parison to current studies on bar- halo interactions. The 
best fit ellipsoid to the bar with aj, — 12.5 kpc has axis 
ratio ai : ai : a% = 3.6 : 1.4 : 1 The mass in the disk 
within the elliptical contour of q = 0.4 with bar length 



14 



at, = 12.5 kpc is Mb = 3.7 x 10 M Q compared with a 
total disk mass of M d = 5.5 x 10 10 M Q . The bar there- 
fore represents 2/3 of the total disk mass. The mass of 
the DM halo within the sphere of radius ab = 12.5 kpc 
is Mh(r < a&) = 6.1 x 10 10 M Q compared with a total 
halo mass M h . to t = 3.0 x 10 11 M Q . The ratio of bar 
to total halo mass enclosed is M{,/M^(r < at,) ~ 0.6. 
We will see below that this perturbation has a small ef- 
fect on the density profile of the DM halo and is not 
sufficient to create a flat density core as seen in recent 
work with a rigid bar evolving in a spherical A^-body 
halo with M^/Mh = .5 — 1.0 and a thinner b ar with 
a x : a 2 : a 3 = 10 : 2 : 1 (IWeinberg fe Katzll2007bh . 

4. HALO DENSITY PROFILE EVOLUTION 

Our next task is t o examine the evolution o f the dark 
halo density profile. IWeinberg fc Kat d ((2002) originally 
demonstrated that a thin rigid bar rotating within a 
cuspy dark halo can disturb the central density pro- 
file and set up a constant density core and follow-up 
work with improved methods and resolution confirmed 
that result for their parti c ular choice of ba r parameters 
(|Weinberg fc Katdl2007blh ISellwoodl (|200ll has recently 
verified these results using independent methods but has 
questioned the applicability of the results of the dynam- 
ics of an idealized thin, rigid bar to real barred galaxies. 
The model described here differs from these studies by 
examining a self- consistent model of a bar forming from 
an instability in an exponential disk within a cusped dark 
halo and so arguably represents a system closer to real- 
ity. A detailed characterization of the model here will 
allow us to compare our results to these other studies. 

Figure [21] presents the evolution of the density profile 
as a function of mass resolution. The plots show the pro- 
files at changing times along with the differential change 
with respect to the initial profile. Gravitational softening 
introduces an artificial density core within ~ 3 Plummcr 
softening lengths but beyond this radius the plots clearly 
show similar behavior in the density profiles. A compar- 
ison of the final density profiles at different resolutions 
again shows similar behavior beyond the softening ra- 
dius and convergence to a similar central behavior. The 
central density profile actually increases by 1.7 x while 
maintaining a central cusp (Fig. |2"2"|) . The likely cause 
of this increase is the halo respon se to the forming bar 
(|Sellwoodl 12001 IColm et al.ll2006h . Once the bar buck- 
les it forms a more concentrated mass distribution in the 
center of the disk and the halo responds by contracting 
adiabatically. In the multi-mass halo with Nh = 10 8 , 
the density cusp is present down to r k 100 pc where 
gravitational softening effects start to influence the dy- 
namics. Within this radius, the halo is well-sampled by 
more than 6000 particles and flattens out into a constant 
d ensity core domina t ed by s oftened gravity. 

IWeinberg fc Katzl (|2007bft (WK herein) have recently 
shown that massive bars can decrease the central density 
of DM halos and disrupt the cusp over a Hubble time in 
some cases at radii of about 20% of the bar length. Our 
bar has a length ab = 12.5 kpc so we should expect to 
see distortions of the density profile at r « 2 kpc while 
in fact we see no signs of a density core developing until 
softened gravity dominates at r = 0.1 kpc in our highest 
resolution case. The results presented here seem to be in 
contradiction so what's going on? The reasons for dis- 



agreement can be understood by comparing the detailed 
properties of the bars used in their models to our self- 
consistently evolved A^-body bar. The WK models are 
rigid, homogeneous ellipsoids of various masses, lengths 
and axis ratios rotating within a live, isotropic A^-body 
halo. Their fiducial model which strongly modifies the 
halo inner profile has a bar length equal to the NFW 
halo scale radius r s , i.e., ab/r s — 1.0, a bar mass equal 
to half the halo mass within this radius M^/Mh = 0.5 
and an axis ratio a\ : a,2 '■ as = 10 : 2 : 1. Our halo is 
also NFW-like but not precisely a NFW model due to 
modifications introduced in setting it up with an embed- 
ded disk and changes induced by bar formation. A good 
proxy for r s in our models is the radius r_2 where the 
density power-law slope 7 = d\ogp/d\ogr ss —2 (For an 
NFW model 7 = -2 at r = r s ). Figure [22] shows that 
r_2 ~ 3 kpc initially. At late times, the 7 profile de- 
velops a wiggle so that 7 = — 2 occurs at two different 
radii, but the average of these two radii is r s rs 5 kpc. 
The final bar length is about 12 kpc (Fig. [19|) and so 
<2(,/r_2 = 2.4. The A^-body bar axis ratio measured 
above is a\ : ai : 0,3 — 3.6 : 1.4 : 1 considerably fatter 
than the fiducial model of WK. Finally, the bar-to-halo 
mass ratio within the bar length is Mb/M^ = 0.6. The 
main differences between the WK fiducial model and our 
model is that their halo is more extended and the bar is 
much thinner overall while the mass ratios are compara- 
ble. The WK models with thicker bars with 02/ ai > 0.3 
are the closest ones to our A'-body mode l s and a ccord- 
ing to their Fig. 13 in IWeinberg fc Katzl (|2007bf ) cause 
no appreciable change in the density profile. So we find 
no inconsistency with their most closely matching model. 

While the thin, massive bars described by WK have 
strong effects on halo profiles, the thicker bars that de- 
velop through the recurrent buckling instabilities are 
more relevant to the evolution of real barred galaxies. 
Thin bars are subject to the dynamical buckling in- 
stability and thicken quickly. Moreover, the strongest 
bars, i.e., those with b/a ~ 0.2 show a rapid decrease 
in the phase space available to regular orbits and hence 
an increase in the fraction of chaotic o rbits in the bar 
(|Martinez-Vabuesta fc Shlosmanl I2004D . While verti- 
cally thin rigid bars are immune to any instabilities, the 
DM particle orbits in the cusp can be destabilized by 
the mere presence of a more massive analytical potential 
mixed with the live potential. Vertically thinner bars, 
i.e., smaller c/a, will be more efficient in destabilizing 
the DM trajectories, by analogy with smaller b/a. In 
any case, the c/a = 0.1 thin bars used by WK cannot be 
justified over a Hubble time. They are supported neither 
by observations or high-resolution numerical simulations. 

We conclude that bars that form self-consistently in N- 
body simulations from the instability of an exponential 
disk in NFW-like DM halos do not destroy the density 
cusp and in fact can increase the halo central density 
slightly. Our mass resolution study shows a clear con- 
vergence in behavior to higher resolution and the central 
characteristics of dark halos are limited only by the parti- 
cle softening and diminishing particle numbers. We now 
explore the detailed orbital dynamics of the bars to un- 
derstand angular momentum transport from the bar to 
the halo through low order resonances. 
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Fig. 21. — Evolution of the dark halo density profile for different mass resolutions at t = 0, 2.4, 4.7 and 7.1 Gyr (black, red, green 
and blue lines). The dotted vertical line shows the value of the Plummer softening length for each resolution: e = 50,20,10 and 10 pc for 
Nh = 10 6 , 10 7 , 10 s and 10 8 (multi-mass) halo models respectively. For r > 10 kpc, the density profile does not change significantly. In the 
range 1 < r < 10 kpc, the density increases roughly 1.7x, showing adiabatic contraction in response to the buckling instability and the 
formation of a centrally concentrated bulge-like bar within the disk. The logarithmic slope of the density profile is a ~ — 1 to within a few 
softening lengths from the center. A constant density core develops within the center with a core radius that depends on and e with 
typical values of ~ 5e. As increases and e decreases, central density increases while the core radius declines. The existence of a small 
core is consistent with relaxation due to softened gravity rather than forcing by the bar. 



5. BAR ORBITAL DYNAMICS 

Angular momentum is transferred from the bar 
to the halo through l ow o rd er orbital resonances 
(iLvnden-Bell fc Kalnaii Il972t iTremaine fe Weinberg! 
19841: Weinberg 1985). Following the convention of 
Weinberg fc Katzl ( 2007af ). the condition for planar res- 
onances is l\D. r + ^2^0 = mflb where (li,l2,m) are an 
integer triplet with radial and azimuthal o rbital frequen- 
cies p, r and 0^, and bar patter n speed (jAthanassoulal 
l2002HWeinberg fc Kat zll2007ar i. In the discussion below, 
we also use the parameters = and k = Q r to refer 
to the true orbital frequencies rather than the epicyclic 
approximations. Bars are predominantly a m = 2 distur- 
bance so integer pairs l\ : I2 with m — 2 correspond to 



various resonances with the more important ones being 
the inner and outer Lindblad resonances (ILR —1:2 and 
OLR 1 : 2) and the corotation resonance (COR 0:2). 
Other important resonances that may transfer angular 
momentum occur with I2 = —2,0 including the direct 
radial resonance (DRR 1:0) discussed by WK. 

We focus our analysis on resonances with I2 = 2 that 
are responsible for the bulk of angular momentum trans- 
fer. A simple way of characterizing the low order res- 
onances is with the d i mensi o nless frequency 77 = (Q — 
p,h)/K (jAthanassoulal 120021 : iMartinez-Valpuesta et al.l 
|2006| ). The half integer values of 77 correspond to low 
order resonances with r\ = —1/2, 0, 1/2 corresponding to 
the OLR, COR, and ILR respectively. Most angular mo- 
mentum is transferred to and from orbits that satisfy this 
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log r/kpc 

Fig. 22. — Change in the density power-law index profile 7 = 
(ilogp/dlogr at t =0 Gyr (solid), 4.7 Gyr (dotted) and 7.1 Gyr 
(dashed). The density profile maintains a cuspcd profile with 7 < 
— 1 down to r = 0.1 kpc well within the the scale radius, r_2 ~ 5 
kpc at late times. A constant density core does not develop in 
response to the bar and the halo maintains its cusp to the limit of 
gravitational softening. 

resonant condition. As the bar loses angular momentum 
and £!{, declines, the population of halo particles in res- 
onance with the bar changes. The potential of the halo 
also readjusts in respons e to the bar, so orbital fre quen- 
cies can change as well. IWeinberg &: Kat z (2007a) have 
argued that the resonances may only occur over a small 
fraction of the halo mass so that poorly resolved halos 
may not have sufficient numbers of particles to absorb 
angular momentum. Furthermore, noise in lower resolu- 
tion simulations can cause particles to move in and out 
of resonance in a diffusive manner leading to an incorrect 
determination of angular momentum transfer. They es- 
timate that as many as 10 s halo particles are necessary 
to both populate resonances and suppress noise to con- 
verge on the correct behavior. We examine these effects 
directly at different resolutions by studying the behavior 
of angular momentum transfer and orbital resonances us- 
ing our models at the recommended mass resolution and 
see if the results do converge. 

5.1. Net Angular Momentum Transfer 

We first examine the net angular momentum transfer 
evolution as a function of mass resolution (Fig. [24]) . We 
have again offset the times at different resolution so that 
they are synchronized with the time of maximum | ^4.2 1 - 
The initial behavior is similar though there is no clear 
trend in behavior between resolutions from t = 2.4 — 7.2 
Gyr reflecting the variance from different random initial 
conditions. At late times, however, the rate of angular 
momentum transfer from the bar to the halo depends on 
resolution, with lower resolution simulations transferring 
J more quickly than the highest resolution case. From 
t = 7.2 — 9.4 Gyr the rate of change J is about two 
times larger for N h = IM than N h = 100M. This effect 
could be the result of noise broadening the resonant inter- 
action though this interpretation is complicated by the 
variance in behavior due to different initial conditions. 



In summary, there is a measurable difference in angu- 
lar momentum transfer between high and low resolution 
with the lowest resolution model transferring about 10% 
more angular momentum. 

5.2. Halo Orbital Resonances 

We quantify the importance of low order resonances 
for angular momentum transfer in our models us- 
ing a modi fied version of the orbi tal spectral analy- 
sis method (|Binnev fc Sper gcl 1982) in a frozen r otat- 
ing potential as described by lAthanassoulal (|2002| ) and 
iMartinez-Valpuesta et al.l (|2006l ). We determine the 
principal orbital frequencies k and for a set of N^h 
randomly chosen particles in the halo with Nnh ~ 10 6 for 
simulations with Nh > 10 6 and N^h = Nh for smaller 
simulations. We then compute the potential and force 
field on a grid with variable spacing to be used for in- 
terpolating forces for test particle integrations. Orbital 
frequencies are determined from test particle integra- 
tion of particles orbits in the frozen potential in a ro- 
tating frame at the bar's pattern speed for the time of a 
given snapshot. The orbits are integrated for about 50 
bar rotations, starting at three representative times - 
t = 2.4, 4.7 and 7.1 Gyr for the N h = 10 s single and multi 
mass simulations. Appropriate time offsets are applied 
as discussed above to lower Nh simulations to synchro- 
nize the time of maximum IA2I. Each orbit was sampled 
with at least 200 constant timesteps per azimuthal pe- 
riod, and overall by 10K timesteps. The use of constant 
timesteps simplifies Fourier decomposition of the orbit 
time series. Most decompositions lead to a line spec- 
trum allowing easy identification of frequencies Q and K 
though occasionally the spectrum is more complex and 
no frequencies can be uniquely identified. 

We present the results of the spectral analysis for DM 
halo orbits in Fig. [25] for various resolutions at t = 7.1 
Gyr for the 100M particle runs. Again, we account for 
the time offsets discussed above for the lower resolution 
runs for a fair comparison. The particles are binned in 
frequency rj with a bin width An = 0.005. Figure [25] 
shows the distribution of the particle number fraction 
(or mass fraction in the multi-mass model case) as a 
function of the dimcnsionless frequency. The main reso- 
nances - ILR, COR, OLR - are present along with higher 
order ones with COR being the most populated reso- 
nance. The relative height of the peaks begins to con- 
verge when N > 10 6 and the behavior is quite similar. 
The peak bins contain a few percent of the total parti- 
cle numbers or ~ 10 6 particles in the largest case and so 
provide good coverage of the resonance for angular mo- 
mentum transfer. We can define the amount of mass in 
"resonance" as the sum over particles with dimensionless 
frequencies in the range 8n ± 0.05 at half integer values 
of n. When measured this way about 7% of the total 
halo mass is in resonance instantaneously at late times 
when the bar has reached quasi-equilibr ium and is slow- 
ing down. iTremaine fc Weinberg] (|1984l ) speculated that 
orbits may have become trapped in resonance if the bar 
slowed down gradually and we can check whether this 
trapping is significant. A comparison of the particles in 
resonant peaks at t = 4.7 Gyr with those at t = 7.1 Gyr 
shows that only a small fraction migrate between reso- 
nances as the system evolves. Of the 7% of the total 
mass in resonance at t = 4.7 Gyr only 1.5% are still in 
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log r/kpc 

Fig. 23. — A comparison of density profiles at t = 7.1 Gyr for different halo particle numbers N/ x . We also show the initial density profile 
(dashed line) and the best fit NFW model curve (dotted line) to the initial profile over the range < r < fOO kpc. The NFW parameters 
for the fit are r a = 4.3 kpc, v max = f60 km s _1 , where v max is the maximum circular velocity at r = 2.f6r s = 9.3 kpc. Note that this halo 
is more concentrated than the typical galactic dark matter halos in cosmological simulations. We use the NFW formula to characterize the 
profile and show that an r _1 cusp extends to within at least fOO pc of the center. The dotted vertical lines show the softening length e 
used at different resolutions. As Nf, increases, the central density increases and the core radius decreases suggesting that the core behavior 
is due to mass resolution rather than forcing by the bar. 



resonances at t = 7.1 Gyr with most particles moving 
out of resonance. As the bar is braking, new orbits are 
brought into resonance while orbits that have acquired 
angular momentum move out of resonance. In this sense, 
the resonance is broad and a significant fraction of halo 
orbits participate in angular momentum exchange with 
the bar. 

Figures [26] and [27] show the resonant transfer of angu- 
lar momentum between the two snapshots at t = 4.7 Gyr 
and t = 7.1 Gyr. We plot the distribution of the change 
in z-angular momentum AJ Z versus rj measured for the 
particles at t — 4.7 Gyr. Most angular momentum is 
absorbed in the halo at the COR and OLR with smaller 
amounts absorbed at higher order resonances. However, 
some J z is emitted and lost from the ILR in accord with 
fundamental ideas of angular momentum transport in 



stellar systems (]Lvnden-Bell fc Kalnaisl [T972). The dis- 
tributions when viewed with an expanded vertical scale 
show nice convergence in detailed behavior at higher res- 
olution (Fig. EH]). For N h > 10 7 , w 50% of the total 
transferred angular momentum is is in the resonant peaks 
(within Srj = ±0.05 while for Nh < 10 6 we find less than 
30% in the peaks with this same definition. The lower 
resolution simulations are clearly more susceptible to dif- 
fusion. Nevertheless, despite these differences the total 
angular momentum transferred is similar for N > 10 6 
suggesting that the diffusive process that broadens reso- 
nances is not a serious problem for the global evolution 
of the system. 
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Fig. 24. — Evolution of the net angular momentum in the disk 
and halo at different mass resolution. Total angular momentum is 
conserved to within 1%. The evolution is similar for all resolutions 
plotted starting at Nf l = 10 6 . The rate of angular momentum 
transfer is slightly smaller at later times for higher resolution sim- 
ulations leading to about a 10% difference in the total amount of 
transferred angular momentum in the multi-mass case suggesting 
a significant but small affect due to resolution. 
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Fig. 26. — Net change in the DM halo particle angular momen- 
tum between t = 4.7 and 7.1 Gyr for particles binned as a function 
of the dimensionless frequency rj measured at t = 4.7 Gyr. The ma- 
jority of angular momentum is gained through the CR resonance 
at rj = though some angular momentum is lost at the ILR at 
rj = 0.5. The peaks are sharper at higher resolution. 
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Fig. 25. — Distribution of DM halo particles as a function of 
the dimensionless frequency rj. Resonant spikes at the half integer 
values of rj correspond to low order resonances. The bin width is 
Sij = 0.005. The distributions are similar as a function of mass 
resolution. 



5.3. Resonances in Phase Space 

Finally, we examine the change in halo phase space 
density by computing the particle number density in 
(E, J z ) space and computing the difference between t = 
and t = 150 in model mlOOM in a similar way to 
Holl ev-Bockelmann et~aT1 (12005). In this way, we clearly 
see the resonant regions visible as discrete islands of 
particle overdensity in (E,J Z ) space (Fig. 125)1 . We can 



Fig. 27. — Same as Fig. 1261 with the vertical scale expanded by 
10 X. The detailed distributions of the change in angular momen- 
tum are similar between the peaks at higher resolution. 



also overplot the values of (E, J z ) for the particles found 
in the resonant spikes in the analysis at the final time 
t = 150 to see where they lie in phase space. Figure [29] 
clearly shows that the peaks in phase space density are 
directly related to the discrete resonances extracted from 
our spectral analysis. An accompanying animation to 
FiglSSlpresents the time evolution of the differential num- 
ber density in phase space and reveals how the resonant 
islands move through a large fraction of the halo mass. 
By counting particles in resonant peaks at different times 
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Change in particle number density in (E, J z ) space 
and t = 150 (7.0 Gyr) for the N h = 10 8 single 
The resonant regions show up clearly as peaks (red 
regions) in phase space in the left panel. The blue-black region is 
a valley where a halo bar rotating along with the disk bar and so 
de-populated the negative J z of phase space at the ILR. See Video 
4 to view the time evolution of the particle phase-space density. 
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Fig. 29. — On the phase number density map, we overplot 
the (E, J z ) coordinates of a subset of particles located at discrete 
resonances at t = 150 within 8r] = ±0.05 (black-ILR-r; = 0.5, 
red-COR-?; = 0.0, green-OLR-jj = —0.5, blue-?} = —1.0, magenta- 
r\ = —1.5, and cyan-7/ = —2.0. The resonant particles lay directly 
on top of the peaks and so identify the specific resonant regions in 
phase space. 



we estimate that roughly 20-30% of the halo particles are 
in resonance with the bar at some time in their history. 
Since such a large fraction of particles are involved in 
angular momentum transfer then even lower resolution 
simulations can do a reasonable job of following the evo- 
lution of the bar. 

6. CONCLUSIONS 

We have carried out a comprehensive set of experi- 
ments to explore the evolution of a self-consistent bar in 
a galactic model with an exponential disk and cuspy DM 
halo using resolutions with 10 4 ~ 8 DM particles and a 
single experiment using a multi-mass method with an ef- 
fective resolution of 10 10 . Our highest resolution exceeds 
by far the level prescribed by IWeinberg fc Kat d (|2007a| ) 
necessary to achieve convergent behavior in bar galaxy 
dynamics. We have applied various diagnostics of bar 
evolution as a function of mass resolution including bar 
growth, pattern speed evolution, halo density cusp evo- 
lution and the resonant transfer of angular momentum 
from the bar to the DM halo. In almost all cases, the gen- 
eral behavior is similar at most but the lowest resolutions 
with the convergence occurring around 10 6-7 , depending 
on the phenomenon. Sellwood (2008) has also explored 
similar effects in a mass resolution study with rigid bars 
in cuspy spherical halos with ~ 10 8 particles and come 
to similar conclusions about minimal resolution require- 
ments. Notably, in this model the density cusp is not 
destroyed by the formation of the bar in apparent con- 
tradiction to the results of WK. Our best explanation is 
that the thick bar that form in our self-consistent models 
has a weaker affect than the rigid thin bars in the work 
of WK and we question the applicability of these thin 
bar models over a Hubble time in light of the buckling 
instability. 



The strongest argument for convergence comes from 
the spectral analysis of orbits in the rotating barred po- 
tential at different resolutions that shows in detail similar 
distributions as a function of the dimensionless frequency 
rj both in mass fractions and angular momentum trans- 
ferred between different times. Analysis of the change 
in phase space density show that resonant islands sweep 
through the phase space as the bar loses angular mo- 
mentum leading to effectively broader resonances with 
as much as 20-30% of the halo mass absorbing angular 
momentum from the bar. 

Future studies should examine the bar instability self- 
consistently using the same initial conditions with dif- 
ferent iV-body methods to resolve current inconsistent 
results on the cusp/core evolution of DM halos as well 
as explore detailed behavior in phase space. The model 
snapshots and initial conditions from this study are freely 
available to researchers in the area who wish to verify our 
results against their own codes and methods. 
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